################################################################################
##### R Code To Replicate "Honor and War"        ###############################
##### Authors: Allan Dafoe and Devin Caughey     ###############################
##### NB: This takes several days to run in full ###############################
################################################################################

################################################################################
#### SET-UP ####################################################################
################################################################################

### Clear workspace
rm(list=ls())
gc()

### Define home directory
setwd(home.dir <- "~/Dropbox/Shared/Work/Honor/Replication-Final")
##setwd(home.dir <- "/nfs/projects_nobackup/o/oldpolls/Honor/Replication-Final")

### Load libraries
stopifnot(require(foreign))
stopifnot(require(Matching)) 
stopifnot(require(rgenoud))
stopifnot(require(lme4)) 
stopifnot(require(coxme)) 
stopifnot(require(sfsmisc))
stopifnot(require(permute)) 
stopifnot(require(plyr)) 
stopifnot(require(dplyr))
stopifnot(require(coin)) 
stopifnot(require(matlab)) 
stopifnot(require(ggplot2))
stopifnot(require(ordinal)) 
stopifnot(require(RItools))
stopifnot(require(reshape2)) 
stopifnot(require(car))
stopifnot(require(TeachingDemos)) 
stopifnot(require(NPC))

ToDate <- function (numeric.date) {
  as.Date(numeric.date, origin = "1970-01-01")
}
mysw <- function (...) {
  setwd(paste(..., sep = "/"))
}
FileName <- function (path=NULL, name=NULL, ext=NULL, replace=FALSE) {
  p <- paste(path, collapse = "")
  n <- paste(name, collapse = "")
  e <- paste(".", ext, sep = "")
  d <- format(Sys.Date(), "%y%m%d") 
  for (l in seq_along(letters)) {
    if (l > 1) old_file_name <- file_name
    file_name <- paste0(p, d, paste(n, letters[l], sep="-"), e)
    if (!any(grepl(file_name, list.files()))) {
      if (replace) file_name <- ifelse(l > 1, old_file_name, file_name)
      break
    }
  }
  cat("\nFile name:", file_name, "\n\n")
  return(file_name)
}
SourceLatest <- function (name, dir=".", ...) {
  all.files <- list.files(dir)
  files <- all.files[grepl(name, all.files) & grepl("\\.R\\>", all.files)]
  dates <- sub(".*([0-9]{6}).*", "\\1", files)
  file <- paste(dir, files[which.max(dates)], sep="/")
  cat("\nOpening:", files[which.max(dates)], "\n")
  source(file, ...)
}

mysw(home.dir, "4-Analysis")
##SourceLatest("Analysis", echo = TRUE, max.deparse.length = Inf)

mysw(home.dir, "4-Analysis/out")
txt.out <- FileName("Honor",, "out")
txtStart(file=txt.out)
getwd()

### Load data (four South codings)
mysw(home.dir, "3-Analysis-Data")
## 0: G. W. Bush Southern, Truman non-Southern (original coding)
mid1 <- read.dta("15-03-13-USMIDs1.dta")
## 1: G. W. Bush Southern, Truman Southern
mid2 <- read.dta("15-03-13-USMIDs2.dta")
## 2: G. W. Bush non-Southern
mid3 <- read.dta("15-03-13-USMIDs3.dta")
## 3: G. W. Bush non-Southern, Truman Southern
mid4 <- read.dta("15-03-13-USMIDs4.dta")

## For storing results
match.ls <- list()
balance.ls <- list()
result.df <- data.frame(SouthCoding = character(1),
                        MID_subset = character(1),
                        cov_set = character(1),
                        replace = character(1),
                        test_stat = character(1),
                        ForceUS_T0 = numeric(1),
                        DurationUS_T0 = numeric(1),
                        OutcomeUS_T0 = numeric(1),
                        ForceUS_p = numeric(1),
                        DurationUS_p = numeric(1),
                        OutcomeUS_p = numeric(1),
                        ForceUS_p_adj = numeric(1),
                        DurationUS_p_adj = numeric(1),
                        OutcomeUS_p_adj = numeric(1),
                        global_p = numeric(1))

mysw(home.dir, "4-Analysis/save")
result.file <- FileName("Results", "", "dta")

################################################################################
##### Loop over South codings ##################################################
################################################################################
SouthCode <- 1
for (SouthCode in SouthCode:4) {
  SouthLabel <-
    c("original", "Truman in", "Bush out", "Bush out and Truman in")[SouthCode]
  cat("\n*** South Coding:", SouthLabel, "***\n")
  ## MID data for this South coding
  mid <- list(mid1, mid2, mid3, mid4)[[SouthCode]]
  ls()

  ## better labels
  mid$LeaderName <- Recode(mid$Leader, "
'AdamsJ' = 'Adams, J.';
'AdamsJQ' = 'Adams, J. Q.';
'BushGHW' = 'Bush, G. H. W.';
'BushGW' = 'Bush, G. W.';
'HarrisonB' = 'Harrison, B.';
'HarrisonWH' = 'Harrison, W.';
'JohnsonA' = 'Johnson, A.';
'JohnsonLB' = 'Johnson, L.';
'RooseveltFD' = 'Roosevelt, F.';
'RooseveltT' = 'Roosevelt, T.';
'Van_Buren' = 'Van Buren'
")

### Define subsets of data
  mid <- mutate(mid,
                EndedWithinTerm = Pres2 == "",
                ## Non-fishing MIDs
                nonfishing = FishDisp == 0 & !is.na(FishDisp),
                ## Non-fishing MIDs where both countries were originators
                primary = Primary == 1 & FishDisp == 0 & !is.na(Primary),
                ## Non-fishing MIDs where both countries were alone on their side
                bilateral = BothAlone == 1 & FishDisp == 0 & !is.na(BothAlone),
                ## Presidents covered by time range of MID dataset (1816-2010)
                in.data = TermEndDate1 > as.Date("1816-01-01") &
                    TermStartDate1 < as.Date("2010-06-01"),
                ## Presidents who experienced at least one non-fishing MID
                non0nofish = nMIDsNoF > 0 & !is.na(nMIDsNoF),
                ## Presidents who experienced at least one primary MID
                non0prim = nMIDsPrim > 0 & !is.na(nMIDsPrim),
                ## Presidents that experienced at least one bilateral MID
                non0bilat = nMIDsAlone > 0 & !is.na(nMIDsAlone),
                ## First presidential observation in each MID subset
                lf = LeaderOb == 1 & non0nofish,
                lp = LeaderOb == 1 & non0prim,
                lb = LeaderOb == 1 & non0bilat,
                SouthFactor = factor(ifelse(South1 == 0, 'Non-Southern',
                    'Southern')))

################################################################################
##### CREATE VARIABLES #########################################################
################################################################################

  mysw(home.dir, "4-Analysis")
  SourceLatest("CreateVars", echo = TRUE, max.deparse.length = Inf)

################################################################################
#### DATA FRAMES OF VARIABLE SETS ##############################################
################################################################################

  ## Response variables
  Y.df <- subset(mid,, c(ForceUS, DurationUS, OutcomeUS))

  ## Placebo outcomes
  placebo5 <- subset(mid,, c(ForceUS.prop.prev5, ave.duration.prev5,
                            ave.outcome.prev5))
  placebo10 <- subset(mid,, c(ForceUS.prop.prev10, ave.duration.prev10,
                             ave.outcome.prev10))
  placebo5log <- subset(mid,, c(ForceUS.prop.prev5, ave.log.duration.prev5,
                               ave.outcome.prev5))
  placebo10log <- subset(mid,, c(ForceUS.prop.prev10, ave.log.duration.prev10,
                                ave.outcome.prev10))
  log.placebo <- data.frame(placebo5log, placebo10log)

### Covariate data frames
  ## `basic.covs' (including term length)
  basic.covs <- subset(mid,, c(great.power, super.power, TermStartYear,
                              ongoing.mid, log.yrs.since.war, nMIDs.prev5,
                              TermYrsAdj))
  ## `struct.covs' for structural covariates (including placebos)
  struct.covs <- data.frame(basic.covs, log.placebo,
                            subset(mid,, c(ongoing.war, log.war.dead.pct,
                                          elitevet.imput, South1.prev,
                                          nMIDs.prev10)))
  ## `party.covs' (no Repub bc collinear and few Southerners)
  party.covs <- data.frame(struct.covs, subset(mid,, c(DemDR, Whig)))
  ## recession, unified gov't
  inTerm.covs <- data.frame(party.covs, subset(mid,, c(proprec, unified)))
  ## military
  full.covs <- data.frame(inTerm.covs,
                          subset(mid,, c(military.experience, nMIDsPrim)))

################################################################################
##### PLOT RAW DATA ############################################################
################################################################################

  set.seed(1)
  mid <- mutate(mid, jit = rnorm(nrow(mid), 0, 0.1),
                TermMidDate = as.Date((as.numeric(TermStartDate1)
                    + as.numeric(TermEndDate1)) / 2, origin = "1970-01-01"))

  (mid.plot <- ggplot(subset(mid, nonfishing))
   + geom_rect(aes(xmin = TermStartDate1, xmax = TermEndDate1,
                   ymin = 0.5, ymax = 5.5, fill = SouthFactor))
   + geom_segment(aes(x = TermStartDate1, xend = TermStartDate1,
                      y = 0.5, yend = 5.5), size = .1)
   + geom_segment(aes(x = TermEndDate1, xend = TermEndDate1,
                      y = 0.5, yend = 5.5), size = .1)
   + geom_text(aes(x = TermMidDate, y = 5.52, label = gsub('_', ' ', LeaderName)),
               angle = 90, hjust = 0, size = 2, data = subset(mid, lb))
   + scale_fill_manual(values = c('white', 'grey80'))
   + scale_x_date(breaks = as.Date(c('1816-01-01', '1896-01-01', '1945-01-01',
                      '2011-01-01')), labels = scales:::date_format('%Y'))
   + scale_y_continuous(breaks = 1:5,
                        labels = c('None', 'Threat', 'Show', 'Use', 'War'))
   + scale_shape_manual(values = c(0, 1, 2))
   + scale_alpha_continuous(range = c(0.5, 1))
   + labs(x = NULL, y = 'US Hostility Level', fill = 'President')
   + theme_bw()
   + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
   + expand_limits(y = 6)
   + guides(alpha = FALSE)
   )

  mysw(home.dir, "5-Figures")
  pdf(file = FileName('NofishMIDs1-', '', 'pdf', replace = TRUE),
      width = 9, height = 5)
  print(mid.plot
   + geom_point(aes(x = StartDateUS, y = hostlev + jit, shape = outcome,
                    color = outcome, alpha = BothAlone), size = 3)
   + geom_point(aes(x = StartDateUS, y = hostlev + jit, alpha = BothAlone),
                size = 0.5)
   )
  dev.off()

  mysw(home.dir, "5-Figures")
  pdf(file = FileName('NofishMIDs2-', '', 'pdf', replace = TRUE),
      width = 9, height = 5)
  print(mid.plot
   + geom_segment(aes(y = hostlev + jit, yend = hostlev + jit, alpha = Primary,
                      x = StartDateUS, xend = EndDateUS, color = outcome),
                  size = .5, show_guide = FALSE)
   + geom_point(aes(x = EndDateUS, y = hostlev + jit, alpha = Primary,
                    shape = outcome, color = outcome), size = 3)
   + geom_point(aes(x = StartDateUS, y = hostlev + jit, alpha = Primary),
                size = 0.5)
   ) 
  dev.off()

  print(pres.prim.df <- group_by(subset(mid, Primary == 1), LeaderName) %>%
       summarise(Region = gsub("ern", "", unique(as.character(SouthFactor))),
                 MIDs = n(),
                 Force = mean(ForceUS),
                 Duration = mean(DurationUS),
                 Outcome = mean(OutcomeUS))
   )
  pres.prim.melt <- melt(pres.prim.df, id.vars = c("LeaderName", "Region", "MIDs"))

  mysw(home.dir, "5-Figures")
  pdf(file = FileName('PresSummaryPrimary', '', 'pdf', replace = TRUE),
      width = 9, height = 5)
  print(ggplot(pres.prim.melt)
   + facet_wrap(~variable, scales = "free_y")
   + aes(x = Region, y = value)
   ## + geom_text(aes(label = label, alpha = alpha), hjust = 0, size = 2)
   + geom_boxplot(##aes(weight = MIDs),
       outlier.colour = "white")
   + geom_point(aes(size = MIDs), shape = 1, alpha = 1/3)
   ## + stat_summary(aes(y = value, w = MIDs), fun.y = "weighted.mean",
   ##                geom = "point", size = 3, shape = 5)
   + stat_summary(aes(y = value, w = MIDs), fun.y = mean,
                  geom = "point", size = 4, shape = 5, color = "black")
   + scale_size(range = c(2, 10))
   + theme_bw()
   + theme(panel.grid.major = element_blank(),
           panel.grid.minor = element_blank())
   + labs(x = NULL, y = NULL, size = "# MIDs")
   )
  dev.off()

  (pres.alone.df <- group_by(subset(mid, BothAlone == 1), LeaderName) %>%
       summarise(Region = unique(as.character(SouthFactor)),
                 MIDs = n(),
                 Force = mean(ForceUS),
                 Duration = mean(DurationUS),
                 Outcome = mean(OutcomeUS))
   )
  pres.alone.melt <- melt(pres.alone.df, id.vars = c("LeaderName", "Region", "MIDs"))
  pres.alone.melt <- ddply(pres.alone.melt, ~Region + variable + value, mutate,
                           label = ifelse(length(LeaderName) > 5,
                               paste(as.numeric(LeaderName), collapse = "\n"),
                               paste(as.numeric(LeaderName), collapse = ", ")),
                           alpha = 1/n())
  mysw(home.dir, "5-Figures")
  pdf(file = FileName('PresSummaryBilateral', '', 'pdf', replace = TRUE),
      width = 9, height = 5)
  print(ggplot(pres.alone.melt)
   + facet_wrap(~variable, scales = "free_y")
   + aes(x = Region, y = value)
   ## + geom_text(aes(label = label, alpha = alpha), hjust = 0, size = 2)
   + geom_boxplot(##aes(weight = MIDs),
       outlier.colour = "white")
   + geom_point(aes(size = MIDs), shape = 1, alpha = 1/3)
   ## + stat_summary(aes(y = value, w = MIDs), fun.y = "weighted.mean",
   ##                geom = "point", size = 3, shape = 5)
   + stat_summary(aes(y = value, w = MIDs), fun.y = mean,
                  geom = "point", size = 4, shape = 5, color = "black")
   + scale_size(range = c(2, 10))
   + theme_bw()
   + theme(panel.grid.major = element_blank(),
           panel.grid.minor = element_blank())
   + labs(x = NULL, y = NULL, size = "# MIDs")
   )
  dev.off()


  set.seed(1)
  (ggplot(data = mutate(pres.alone.df, Outcome = jitter(Outcome, 20)))
   + geom_text(aes(x = Force, y = Outcome, size = MIDs^(1/5), label = LeaderName,
                   color = Region), alpha = 1)
   + scale_size(range = c(2, 6))
   + theme_bw()
   + scale_x_continuous(breaks = c(0, .5, 1), limits = c(-.1, 1.1))
   + theme(panel.grid.major = element_blank(),
           panel.grid.minor = element_blank())
   + guides(size = FALSE, color = FALSE)
   )

  set.seed(1)
  (ggplot(data = mutate(pres.alone.df, Outcome = jitter(Outcome, 20)))
   + geom_text(aes(x = Duration, y = Outcome, size = MIDs^(1/5),
                   label = LeaderName, color = Region), alpha = 1)
   + scale_size(range = c(2, 6))
   + theme_bw()
   + scale_x_continuous(breaks = c(0, 100, 200, 300), limits = c(-20, 300))
   + theme(panel.grid.major = element_blank(),
           panel.grid.minor = element_blank())
   + guides(size = FALSE, color = FALSE)
   )

################################################################################
##### ALTERNATIVE EXPLANATIONS #################################################
################################################################################

### Are Southern presidents more likely to get into disputes with non-white
### opponents (racism as alternative explanation)?
  ## Codes for non-white opponents
  nonwhite.countries <- c("COL", "MEX", "UKG", "AUH", "CHL", "JPN",
                          "BRA", "HAI", "CHN", "GUA", "DOM", "MOR",
                          "NIC", "CUB", "TUR", "PAN", "SAL", "HON",
                          "TAW", "PRK", "ECU", "EGY", "SYR", "DRV",
                          "YAR", "CAM", "GUI", "LIB", "IRN", "YPR",
                          "GRN", "IRQ", "SUD", "AFG", "LBR", "VEN",
                          "PER", "PAK", "CDI", "SOM", "INS")
  
  mid$nonwhite <- with(mid, laply(opps, function (op) {
    as.numeric(any(laply(nonwhite.countries, grepl, x = op, fixed = TRUE)))
  }))

  print(summary(lmer(nonwhite ~ South1 + (1 | Leader),
                     data = subset(mid, bilateral))))
  print(summary(lmer(nonwhite ~ South1 + (1 | Leader),
                     data = subset(mid, primary))))
  print(summary(lmer(nonwhite ~ South1 + (1 | Leader),
                     data = subset(mid, nonfishing))))
  
################################################################################
##### PERMUTATION TESTS ########################################################
################################################################################

  HierLogit <- function (y, tr, tl, event, cluster, ...) {
    mod <- glmer(y ~ factor(tr==tl) + (1 | cluster), family = binomial("logit"))
    fixef(mod)[2] / sqrt(vcov(mod)[2, 2])
  }
  HierCoxPH <- function (y, tr, tl, event, cluster, ...) {
    mod <- coxme(Surv(y, event) ~ factor(tr==tl) + (1 | cluster))
    -fixef(mod) / sqrt(vcov(mod)) ## -1 bc coef is for hazard rate
  }
  HierGauss <- function (y, tr, tl, event, cluster, ...) {
    mod <- lmer(y ~ factor(tr==tl) + (1 | cluster))
    fixef(mod)[2] / sqrt(vcov(mod)[2, 2])
  }
  DiffMeanRank <- function(y, tr, tl, ...) {
    Ri <- rank(y)
    return(mean.default(Ri[tr == tl]) - mean.default(Ri[tr != tl]))
  }

  ## test statistics
  mean.tests <- list("DiffMeans", "DiffMeans", "DiffMeans")
  harm.tests <- list("HarmonicWtdMean", "HarmonicWtdMean", "HarmonicWtdMean")
  rank.tests <- list("DiffMeanRank", "LogRank", "DiffMeanRank")
  hier.tests <- list("HierLogit", "HierCoxPH", "HierGauss")
  test.ls <- list(mean = mean.tests, harmonic = harm.tests, rank = rank.tests,
                  hlm = hier.tests)


  ## matching covariate sets
  cov.ls <- list(NULL, log.placebo, basic.covs, struct.covs, party.covs,
                 inTerm.covs, full.covs)
  names(cov.ls) <- c("none", "log.placebo","basic.covs", "struct.covs",
                     "party.covs", "inTerm.covs", "full.covs")

  ## MID subsets
  subset.ls <- list(bilateral = mid$bilateral, primary = mid$primary)

  ## replacement settings
  replace.ls <- list(FALSE, TRUE)

  s <- v <- r <- t <- 1
  for (s in 1:length(subset.ls)) {
    sum(ms <- subset.ls[[s]])
    for (v in 1:length(cov.ls)) {
      summary(mc <- cov.ls[[v]])
      for (r in 1:length(replace.ls)) {
        (rs <- replace.ls[[r]])
        ## UNMATCHED
        if (is.null(mc)) {
          for (t in 1:length(test.ls)) {
            (ts <- test.ls[[t]])        
            cat("\n***", names(subset.ls)[s], names(cov.ls)[v], rs,
                names(test.ls)[t], "***\n")
            if (rs == TRUE || names(test.ls)[t] == "harmonic") next
            npc <- NPC(data = mid[ms, ], test.statistic = ts,
                       alternative = "greater",
                       tr.var = "SouthFactor", tr.label = "Southern",
                       clust.var = "Leader", event.var = "EndedWithinTerm",
                       y.vars = names(Y.df), comb.fun = "NormalCF",
                       FWE.adj = TRUE, na.rm = TRUE, return.matrix = TRUE,
                       n.perms = 1000, seed = 1, print.steps = FALSE)
            print(T0 <- npc$T.matrix["T0", ])
            print(pv <- npc$p.values)
            res <- data.frame(SouthCoding = SouthLabel,
                              MID_subset = names(subset.ls)[s],
                              cov_set = names(cov.ls)[v], replace = NA,
                              test_stat = names(test.ls)[t],
                              ForceUS_T0 = T0[1], DurationUS_T0 = T0[2],
                              OutcomeUS_T0 = T0[3], ForceUS_p = pv[1],
                              DurationUS_p = pv[2], OutcomeUS_p = pv[3],
                              ForceUS_p_adj = pv[4], DurationUS_p_adj = pv[5],
                              OutcomeUS_p_adj = pv[6], global_p = pv[7])
            rownames(res) <- NULL
            mysw(home.dir, "4-Analysis/save")
            if (result.file %in% list.files()) {
              old.results <- read.dta(result.file)
              result.df <- rbind(old.results, res)
            } else {
              result.df <- res
            }
            write.dta(result.df, file = result.file, convert.factors = "string")
            flush.console()
          }
        }
        ## MATCHED
        if (!is.null(mc)) {
          cat("\n***", names(subset.ls)[s], names(cov.ls)[v], rs, "***\n") 
          if (names(subset.ls)[s] == "bilateral") l1 <- mid$lb
          if (names(subset.ls)[s] == "primary") l1 <- mid$lp
          ## Find optimal matches in terms of covariate balance
          set.seed(1)
          gen.out <- GenMatch(Tr = mid$South1[l1], X = mc[l1, ],
                              estimand = "ATT", wait.generations = 10,
                              pop.size = 10000, replace = rs, print = 1, M = 1)
          ## Create match object
          mat.out <- Match(Y = NULL, Tr = mid$South1[l1], X = mc[l1, ],
                           estimand = "ATT", Weight.matrix  =  gen.out, M = 1,
                           replace = rs)
          ## New dataset of matched data
          ind.t <- mat.out$index.treated
          ind.c <- mat.out$index.control
          ind.tc <- c(ind.t, ind.c)
          obs.t <- numeric()
          obs.c <- numeric()
          pair.t <- numeric()
          pair.c <- numeric()
          for (i in seq_along(ind.t)) {
            lti <- mid$Leader == mid$Leader[l1][ind.t][i]
            lci <- mid$Leader == mid$Leader[l1][ind.c][i]
            pair.t <- c(pair.t, rep(i, sum(lti & ms)))
            pair.c <- c(pair.c, rep(i, sum(lci & ms)))
            obs.t <- c(obs.t, which(lti & ms))
            obs.c <- c(obs.c, which(lci & ms))
          }
          pair.mat <- c(pair.t, pair.c)
          obs.mat <- c(obs.t, obs.c)
          N.c <- sum(sapply(mid$Leader[l1][ind.c], function (l) {
            sum(mid$Leader[ms] == l)
          }))
          N.t <- sum(sapply(mid$Leader[l1][ind.t], function (l) {
            sum(mid$Leader[ms] == l)
          }))
          print(data.frame(NameS = mid$LeaderName[l1][ind.t],
                           PrimMIDsS = mid$nMIDsPrim[l1][ind.t],
                           AloneMIDsS = mid$nMIDsAlone[l1][ind.t],
                           NameN = mid$LeaderName[l1][ind.c],
                           PrimMIDsN = mid$nMIDsPrim[l1][ind.c],     
                           AloneMIDsN = mid$nMIDsAlone[l1][ind.c]))
          mid.match <- data.frame(mid[obs.mat, ], Pair = factor(pair.mat),
                                  NorthernMatch = mid$Leader[l1][ind.c][pair.mat])
          pres.match <- data.frame(mid[l1, ][ind.tc, ],
                                   Pair = factor(rep(seq_along(ind.t), 2)))
          match.ls[[length(match.ls) + 1]] <- pres.match
          names(match.ls)[length(match.ls)] <-
            paste(SouthLabel, names(subset.ls)[s], names(cov.ls)[v], rs,
                  names(test.ls)[t])
          ## Balance
          balance.ls[[length(balance.ls) + 1]] <-
            MatchBalance(mid$South1[l1] ~ ., data = mc[l1, ],
                         match.out = mat.out, nboots = 1000) 
          names(balance.ls)[length(balance.ls)] <-
            paste(SouthLabel, names(subset.ls)[s], names(cov.ls)[v], rs,
                  names(test.ls)[t])
          NPC(data = mid.match, tr.var = "SouthFactor", tr.label = "Southern",
              block.var = "Pair", clust.var = "Leader",
              event.var = "EndedWithinTerm", y.vars = rep(names(mc), 2),
              test.statistic = rep(c("HarmonicWtdMean", "KS"),
                  each = ncol(mc)), alternative = "two.sided",
              comb.fun = "ProductCF", FWE.adj = FALSE, n.perms = 1000, seed = 1,
              print.steps = FALSE)

          ## Placebo tests
          MatchBalance(mid$South1[l1] ~ ., data = log.placebo[l1, ],
                       match.out = mat.out, nboots = 1000)
          print(xBalance(South1 ~ ., report=c("all"), strata = pres.match$Pair,
                         data = pres.match[, c("South1", names(log.placebo))]))
          NPC(data = mid.match, tr.var = "SouthFactor", tr.label = "Southern",
              clust.var = "Leader", block.var = "Pair",
              event.var = "EndedWithinTerm", y.vars = rep(names(log.placebo), 2),
              test.statistic = rep(c("DiffMeans", "KS"),
                  each = ncol(log.placebo)), alternative = "two.sided",
              comb.fun = "ProductCF", FWE.adj = FALSE, n.perms = 1000, seed = 1,
              print.steps = FALSE)

          ## Number of MIDs
          print(NPC(data = mid.match, tr.var = "SouthFactor",
                    tr.label = "Southern", clust.var = "Leader",
                    block.var = "Pair", event.var = "EndedWithinTerm",
                    y.vars = "ForceUS", test.statistic = "DiffSumObs",
                    alternative = "two.sided", comb.fun = "ProductCF",
                    FWE.adj = FALSE, n.perms = 1000, seed = 1,
                    print.steps = TRUE))
          
          for (t in 1:length(test.ls)) {
            (ts <- test.ls[[t]])        
            cat("\n***", names(test.ls)[t], "***\n")
            npc <- NPC(data = mid.match,
                       tr.var = "SouthFactor", tr.label = "Southern",
                       block.var = "Pair", clust.var = "Leader",
                       event.var = "EndedWithinTerm", y.vars = names(Y.df), 
                       test.statistic = ts, alternative = "greater",
                       comb.fun = "NormalCF", FWE.adj = TRUE, na.rm = TRUE,
                       return.matrix = TRUE, n.perms = 1000, seed = 1)
            (T0 <- npc$T.matrix["T0", ])
            (pv <- npc$p.values)
            res <- data.frame(SouthCoding = SouthLabel,
                              MID_subset = names(subset.ls)[s],
                              cov_set = names(cov.ls)[v], replace = rs,
                              test_stat = names(test.ls)[t],
                              ForceUS_T0 = T0[1], DurationUS_T0 = T0[2],
                              OutcomeUS_T0 = T0[3], ForceUS_p = pv[1],
                              DurationUS_p = pv[2], OutcomeUS_p = pv[3],
                              ForceUS_p_adj = pv[4], DurationUS_p_adj = pv[5],
                              OutcomeUS_p_adj = pv[6], global_p = pv[7])
            rownames(res) <- NULL
            mysw(home.dir, "4-Analysis/save")
            if (result.file %in% list.files()) {
              old.results <- read.dta(result.file)
              result.df <- rbind(old.results, res)
            } else {
              result.df <- res
            }
            write.dta(result.df, file = result.file, convert.factors = "string")
            flush.console()
          } ## test statistics
        } 
      } ## replacement setting
    } ## Matching variables
  gc()
  } ## MID subset
} ## South coding

mysw(home.dir, "4-Analysis/save")
save.image(file=FileName("Honor", "", "RData"))

